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Abstract 


A new multigrid relaxation scheme, lower-upper 
symmetric successive overrelaxation (LU-SSOR) 
scheme, is developed for the steady-state solution 
of the Euler and Navier-Stokes equations. The 
scheme, which is based on central differences, does 
not require flux splitting for approximate Newton 
iteration. Application to transonic flow shows 
that the new method is efficient and robust. The 
vectorizable LU-SSOR scheme needs only scalar 
diagonal inversions. 

Introduction 


The Reynolds numbers for a large airplane are 
of the order of thirty million. Therefore, laminar 
flow in the boundary layer becomes unstable, result- 
ing in turbulent flow over most of the surface of 
the airplane. However, the computational require- 
ments for the simulation of turbulent flow are 
clearly beyond the reach of current computers. The 
first level of approximation - time averaging the 
rapidly fluctuating components - yields the 
Reynolds-averaged Navier-Stokes equations, which 
require a turbulence model for closure. At the 
present time not much is known about the behavior 
of turbulence in separated regions, and this has 
impeded the development of turbulence models for 
complex three-dimensional flows. Since a univer- 
sally satisfactory turbulence model has yet to be 
found, current turbulence models have to be tailored 
to the particular flow. During the last decade the 
feasibility of solving the Navier-Stokes equations 
has been explored but the methods developed so far 
have been too expensive to use in a routine produc- 
tion mode. 

Recently several implicit schemes have been 
developed successfully in conjunction with a mul- 
tigrid method for steady-state solution of the 
unsteady Euler equations. Although the 
alternating-direction implicit scheme could be 
improved to achieve the expected efficiency of the 
multigrid method in two dimensions, 1 its inherent 
limitations in three dimensions suggest alternative 
approaches. ^ An alternative implicit scheme that 
is stable in any number of space dimensions is based 
on lower-upper (LU) factorization. The LU implicit 
scheme proved to be robust and efficient for high- 
speed flows up to Mach 20 as well as for transonic 
flows. 3-5 it was also shown that a symmetric 
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Gauss-Seidel relaxation method for solving the 
unfactored implicit scheme was a variant of the LU 
implicit scheme. 

The Newton iteration method has been investi- 
gated to solve the steady Euler or Navier-Stokes 
equations. Because of the rapid growth of the 
operation count with the number of mesh cells, the 
system was solved indirectly. Jespersen 6 and 
Hemker and Spekreijse 7 applied the symmetric 
Gauss-Seidel method to the Euler equations, while 
MacCormack° applied the line Gauss-Seidel method 
to the Navier-Stokes equations. In this paper an 
efficient multigrid relaxation scheme is developed 
for approximate Newton iteration. The new lower- 
upper symmetric successive overrelaxation (LU-SSOR) 
scheme requires scalar diagonal inversions while 
the Gauss-Seidel method and the LU implicit scheme 
require block matrix inversions. The use of scalar 
diagonal inversions offers the potential for order- 
of-magnitude speedups when large systems of partial 
differential equations must be solved, for example, 
for hypersonic flows with finite-rate chemistry. It 
is desirable that the matrix be diagonally dominant 
to assure the convergence of a relaxation method. 

The new method based on central differences achieves 
this without the flux splitting that substantially 
increases the computational work per cycle. Numeri- 
cal examples include inviscid and viscous transonic 
airfoi Is. 


Governing Equations 

The Navier-Stokes equations represent gas flow 
in thermodynamic equilibrium. Let t, p, E, T, and 
p be time, density, total energy, temperature, and 
pressure; u and v Cartesian velocity components ; 
and x and y Cartesian coordinates. Then for a 
two-dimensional flow these equations can be written 
as 


3W + aF + 36 _ 3F v + 3G v 
at ax ay ax ay 


where W is the vector of dependent variables, and 
F and G are convective flux vectors: 


W = (p, pu, pv, pE)* 

F = [pu, P u 2 + p, pvu, u ( P E + p)]* 
G = [ P v, p u v, pv 2 + p, v(pE + p)]* 


( 2 ) 


Here * denotes the transpose of a matrix. The 
flux vectors for the viscous terms are 
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Here the viscous stresses are 


r xx = 2pU x - I p(u x + v y ) 


T = u(u + V ) 
xy y x ; 


and 


o 

T = 2uV - 4 u ( U + V ) 

yy y 3 x y' 

where u is the coefficient of viscosity and k 
is the coefficient of thermal conductivity. 

The pressure is obtained from the equation of 
state: 

p = p(y - 1) |e - \ (u 2 + v 2 )J (3) 

where y is the ratio of specific heats. 

Semi-Discrete Finite-Volume Method 


A convenient way to assure a steady-state solu- 
tion independent of the time step is to separate 
the space and time discretization procedures. In 
the semi-discrete finite-volume method one begins 
by applying a semi-discretization in which only the 
spatial derivatives are approximated. Using a 
finite-volume method for space discretization allows 
one to handle arbitrary geometries and helps to 
avoid problems with the metric singularities that 
are usually associated with finite-difference meth- 
ods. Finite-volume methods do not require special 
treatment on a composite grid. 9 The scheme 
reduces to a central difference scheme on a Carte- 
sian grid and is second-order accurate in space pro- 
vided the mesh is smooth enough, and uniform flow is 
an exact solution of the difference equations. The 
Gauss theorem is used to evaluate the viscous flux 
terms A® 


Nonlinear Adaptive Dissipation 


Using a central difference scheme when calcu- 
lating flows with discontinuities typically produces 
flowfield oscillations in the neighborhood of shock 
waves, where the pressure gradients are severe. To 
suppress the tendency for spurious odd and even 
point oscillations and to prevent nonphysical over- 
shoots near shock waves, we augment the scheme by 
artificial dissipative terms. The dissipative term, 
which is constructed so that it is third-order accu- 
rate in smooth regions of the flow, is explicitly 
added to the residual. For the density equation, 
for example, the dissipation d has the form 

d i+l/2,j ~ d i-l/2,j + d i , j+1/2 " d i , j-1/2 
where 


cell area S is equivalent to the inverse of the 
transformation Jacobian determinant. Both coeffi- 
cients include a normalizing factor, S-j+ 1/2 i/At> 
proportional to the length of the cell side) and 

e i+l/2 j iS alS ° made P r °P° rtional t0 the normal- 
ized second difference of the pressure: 




3 i+l,j ~ 2p i,j + Pj-lJ 


Vl.j 


2P< 




i-l> j 


(5) 


in the adjacent cells. The third-order terms pro- 
vide background damping of high-frequency modes. 

The first-order terms, which are needed to control 
oscillations in the neighborhood of shock waves, are 
activated when strong pressure gradients in the flow 
are sensed. The dissipative terms for the other 
equations are constructed from similar formulas. 
Increasing the amount of artificial viscosity 
improves the rate of convergence, although too much 
dissipation can hinder this. However, dissipation 
should be as small as possible so that the accuracy 
of the solution is not degraded. The typical amount 
of the third-order terms is almost negligible when 
compared to the physical viscosity. Schemes con- 
structed along these lines combine the advantages 
of simplicity and economy of computation. They have 
also proved robust in calculations over a wide range 
of Mach numbers. ^ For more accurate capturing of 
oblique shock waves in hypersonic flows, a total 
variation diminishing (TVD) scheme 4 can be used. 


LU-SSOR Scheme 


A prototype implicit scheme for a system of 
nonlinear hyperbolic equations such as the Euler 
equations can be formulated as 

W n+1 = W n - BAt[D x F(W n+1 ) + D y G(W n+1 )] 

- (l- B )At[D x F(W n ) + D y G(W n )] (6) 

where D x and Dy are central difference oper- 
ators that approximate 3/ax and a/ay. Here n 
denotes the time level. An enormous number of com- 
putations must be performed when the scheme is in 
this form because coupled nonlinear equations must 
be solved at each time step. Let the Jacobian 
matrices be 


(7) 


- W n 


A = 


aF 


n 

B ‘ aW 


and let the correction be 

«w = w n+1 


The scheme can be linearized by setting 


(2) 

d i+l/2, j = e i+l/2,j (p i+l,j " p i,j ) 

" e i+l/2,j (p i+2,j " 3p i+l,j + 3p i , j " p i-l,j ) (4) 

Here e ( 2 ) and e(^) denote coefficients of 
the second and fourth differences respectively. The 


F(W" +1 ) = F(W n ) + A6W + 0(ii«Wn 2 ) 

G(W n+1 ) = G(W n ) + B«W + 0( ii«Wii2) 

where 0 is the order of the enclosed terms, and 
dropping terms of the second and higher orders to 
yield 
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[I + BAt(D x A + DyB)]6W + At R = 0 (8) 

where I is the identity matrix and R is the 
residual 

R = D x F(W n ) + DyG(W") 

If a constant b = 1/2, the scheme remains second- 
order accurate in time; for other values of b, the 
time accuracy drops to first order. 

The unfactored implicit scheme (Eq. (8)) 
produces a large block banded matrix that can be 
inverted only by performing a great many computa- 
tions. In addition, a large amount of storage is 
required. If a = 1 the scheme reduces to a Newton 
iteration in the limit At > »: 

(D X A + DyB) <5W + R = 0 (9) 

A diagonally dominant form of Eq. (9) 


A = -^(A 
A" = |(A 
B + = |(B 
B' = |(B 

where 

r A > max 
r g _> max 

Here, X/\ and xg represent eigenvalues of 
Jacobian matrices. 


+ r * 1 ' 


r A I ^ 


+ r B 0 


- r B !) 


( 14 ) 
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(15) 


( D A + + D+A + D B + + D*B )«W + R = 0 
x x y y 

can be written as 

A >u - * WVi.j 


(10) 


Subtract Eq. (12) from Eq. (13) to get 

Kj - A Tj> 5M ij - 

- - < A (j - A 7j>‘"u * Kj - B iiKj 


- V“ij * V"u 


B i,j-l 6W i,j-l 


This may be written as 


(16) 


* - Vu * R n ■ 0 


( 11 ) 


By simulating it with backward and forward 
relaxation sweeps, we obtain the symmetric 
successive overrelaxation (SSOR) method, which can 
be written in two steps as 


(D A + + D B + - A - B")«W = (A + + B + - A" - B~)«W* 
x y 

(17) 


where 


<A (j - A ij>% * a i.i,j‘ u w,j * <»;, - B ij>‘"(i 


«W* = (0/ + D*B + A + + B + ) _1 (- R) (18) 

x y 


* B i,j*i*",,j*i* R ti-° ii2 > 


followed by 
+ 


« A u - A ij>‘“id - A i-i,j‘»i-i,j* A T*i,j‘"i*i,j 

* - <,i-i*Vi-i 

* B T,j*i* w *,j*i * "u - 0 < 13 > 


where D and D are backward difference 

x y 

operators and D and D are forward difference 
x y 

operators. Here, two-point operators are used for 
steady flow calculations. A + , A", B + , and B“ are 
constructed so that the eigenvalues of "+" matrices 
are nonnegative and those of matrices are 
nonpositi ve: 


If we take "+" and matrices as given in 
Eq. (14), 


A + - A' = r A I 
B + - B- = r B I 

Thus, Eq. (17) becomes the LU-SSOR scheme for 
approximate Newton iteration 

(D”A + + D”B + - A" - B")(dV + DyB" 

+ A + + B + )«W = - (r A + r B )(D x F + D y G) (19) 

which can be inverted in two steps. 

For the Navier-Stokes equations, F and G 
are replaced by F - F v and G - G v in Eq. (19). 
That is. 
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A + + D y B + - A - B j ^D*A + D y B + A + + B + j 5W 

= - (r. + r D ) To (F - F ) + D (G - G ) 
v A B' L * v' y v v' 

- 4 2) D X W + 4 4) D I W - 4 2) °5 W + 4 4) D y W ] 

( 20 ) 

where d 2 and D 4 denote the second dif- 
ference and the fourth difference operators 
respecti vely. 

Since one-sided difference schemes are natu- 
rally dissipative, no implicit smoothing is required 
on the left side. Only adaptive dissipation terms 
are explicitly added to the residual on the right 
side. It is interesting to note that the present 
numerical method eliminates the need for block diag- 
onal inversions without using the diagonal ization 
process. This is an especially desirable feature 
for the analysis of hypersonic reacting flows. The 
Llf family of algorithms are fully vectorizable along 
i + j = constant lines on a vector computer. 

Multigrid Method 

The underlying idea of a multigrid method is 
to transfer some of the task of tracking the evolu- 
tion of the system to a sequence of successively 
coarser meshes. 11 This has two advantages. 

First, the computational effort per cycle is reduced 
on a coarser mesh. Second, the use of larger con- 
trol volumes on the coarser grids tracks the evolu- 
tion on a larger scale, with the consequence that 
global equilibrium can be more rapidly attained. 

In general one can conceive of a multigrid method 
using a sequence of independently generated coarser 
meshes which are not associated with each other in 
any structured way. Here attention is restricted 
to the case in which coarser meshes are generated 
by eliminating alternating points in each coordinate 
direction. Thus simple rules can be formulated for 
the transfer of data between grids. The cells of 
the fine mesh can be amalgamated into larger cells 
which form a coarser mesh. Then in each coarse 
mesh cell the conservation laws are represented by 
summing the flux balances of its fine mesh cells; 
consequently, the evolution on the coarse mesh is 
driven by the disequilibrium of the fine mesh equa- 
tions. The multigrid method used here is the cell- 
centered method which was used for the implicit 
schemes 

Results 


The first test case was for inviscid transonic 
flow past the NACA 0012 airfoil at 1.25° angle of 
attack. The freestream Mach number was 0.8. Nonre- 
flecting boundary conditions were used to absorb 
the waves impinging on the far-field boundary.2 
Figure 1 shows the plot of Mach number contours. 
Convergence histories of the LU-SSOR scheme are 
compared with those of the LU implicit scheme of 
Ref. 2. Figure 2 shows the logarithm of the aver- 
age density residuals, and Fig. 3 shows the lift 
histories. Five-level multigrid calculations were 
performed on a 128 by 32 C-mesh without grid 
sequencing. Uniform flow was given as the initial 
condition. As the results show, the convergence 
rate of the LU-SSOR scheme is about 30 percent 
faster than that of the LU impMcit scheme. More- 
over, the computational work per cycle for the 


LU-SSOR scheme is about 30 percent less than that 
for the LU implicit scheme since the LU-SSOR scheme 
does not need block diagonal inversion. This offers 
the potential for order-of-magnitude speedups when 
large systems of partial differential equations 
must be solved. For example, the size of Jacobian 
matrices for scramjet flows can be over 20 by 20 in 
three dimensions if 15 species equations are added. 
It seems that slow convergence of the LU-SSOR 
implicit scheme in Ref. 2 is due to the method of 
local time stepping that was used in Ref. 2. 

The next case was for viscous laminar flow past 
the NACA 0012 airfoil at Mach 0.5, Reynolds number 
5000, and zero angle of attack. The adjabatic wall 
boundary condition was used at the body surface. 
Calculations were performed on a stretched 192 by 
48 C-mesh. A convergence acceleration technique 
such as enthalpy damping was not used for the 
viscous flow calculations. Figure 4 shows the Mach 
number contours while Fig. 5 shows velocity vectors. 
Convergence history with a two-level multigrid is 
shown in Fig. 6. Theoretically, improved conver- 
gence rates might be obtained by using more levels 
of multigrid. 

The last case was for viscous turbulent flow 
past the RAE 2822 airfoil attach 0.73, Reynolds 
number 6.5 million, and 2.79° angle of attack. The 
Reynolds-averaged Navier-Stokes equations were 
solved using a Baldwin-Lomax turbulence model. ^ 
Transition was fixed at 3 percent chord. Mach 
number contours are shown in Fig. 7 (the dashed line 
denotes the sonic line). Convergence histories on a 
192 by 48 mesh with two-level multigrid are shown: 
Fig. 8 shows the residual history and Fig. 9 shows 
the lift and drag history. The convergence based on 
the average residual is continuous although it is 
slowed after dropping about three orders of magni- 
tude. It seems possible to improve the accuracy 
and the convergence rate further by using a better 
distribution of grid points. 

Conclusion 

A new relaxation scheme derived for approxi- 
mate Newton iteration is applied to both Euler and 
Navier-Stokes equations. The LU-SSOR scheme is not 
only stable in three dimensions but also fully vec- 
torizable. Transonic airfoil calculations confirm 
that the scheme is robust and efficient. Moreover, 
the LU-SSOR scheme is promising for high-speed 
reacting flow calculations. 
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Figure 2.- Convergence history for inviscid transonic flow. 
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Figure 3.- Lift history for inviscid transonic flow. 



Figure 5.- Velocity vectors for viscous laminar flow. 



Figure 6.- Convergence history for viscous laminar flow. 
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Figure 8.- Convergence history for viscous turbulent flow. 


LIFT AND DRAG 



Figure 9. - Lift and drag histories for viscous turbulent 
flow. 
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